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Abstract. 

We study the electronic structure of graphene with a single substitutional 
vacancy using a combination of the density-functional, tight-binding, and impurity 
Green's function approaches. Density functional studies are performed with the 
all-electron spin-polarized linear augmented plane wave (LAPW) method. The 
three sp^cr dangling bonds adjacent to the vacancy introduce localized states 
(Vcr) in the mid-gap region, which split due to the crystal field and a Jahn- Teller 
distortion, while the pzTr states introduce a sharp resonance state (Vtt) in the 
band structure. For a planar structure, symmetry strictly forbids hybridization 
between the a and the n states, so that these bands are clearly identifiable in 
the calculated band structure. As for the magnetic moment of the vacancy, the 
Hund's-rule coupling aligns the spins of the four localized Vcti fj,, V(T2 ti ^'nd 
the Vtt f electrons resulting in a S = 1 state, with a magnetic moment of 2/^^, 
which is reduced by about 0.3^b due to the anti- ferromagnetic spin-polarization 
of the n band itinerant states in the vicinity of the vacancy. This results in the 
net magnetic moment of 1.7 ^b- Using the Lippmann-Schwinger equation, we 
reproduce the well-known ^ 1/r decay of the localized Vtt wave function with 
distance and in addition find an interference term coming from the two Dirac 
points, previously unnoticed in the literature. The long-range nature of the Vtt 
wave function is a unique feature of the graphene vacancy and we suggest that this 
may be one of the reasons for the widely varying relaxed structures and magnetic 
moments reported from the supercell band calculations in the literature. 



PACS numbers: 81.05.Uw; 73.22.-f 



Submitted to: New J. Phys. 



Electronic structure of the substitutional vacancy in graphene: Density-functional and Green's function studies2 



3 



-6- 



''7' 



nil 



DOS 



-VOil 



Figure 1. Sketch of the electronic structure for an isolated substitutional vacancy 
in graphene. The continuum tt and rr bands are shown as dashed curves, while the 
vacancy-induced localized states, Vcr and Vtt, are denoted by straight lines. Ep is 
the Fermi energy. The occupied vacancy states are indicated by solid circles with 
a corresponding net magnetic moment of 2iib- The circular density-of-states 
(DOS) in the midgap region, labelled vr; fj-i indicates schematically the anti- 
ferromagnetic spin-polarization of the tt electron states in the local neighborhood 
of the vacancy. This spin polarization is responsible for the reduction of the 
localized magnetic moment from 2fiB {S = 1) to about 1.7 in our density- 
functional calculation. 



1. Introduction 

Graphene is a material of considerable interest on account of its unusual linearly- 
dispersive Dirac band structure and particle- hole symmetry. [U [2] Vacancy constitutes 
an important defect center, the electronic structure of which forms the basic foundation 
for the understanding of the behavior of more complex defects including impurities. 
Recently it has been suggested that transition-metal doped graphene with vacancies 
may have potential application in hydrogen storage. [3 Experimentally, vacancies in 
graphene have been created intentionally by irradiating materials with electrons and 
ions[ll O [7] and they may also occur in small concentration during the growth 
process. [S] While an ideal graphene sheet is non-magnetic, experimental observation 
of magnetism in carbon systems has been long explained in terms of a variety of defects 
including isolated vacancies, vacancy clusters, or presence of internal or external 
boundaries as in nanoribbons.[7. .Qt IIOI fTTj 

There have been several theoretical studies of the isolated vacancy in graphene 
from first-principles density-functional theory (DFT) [T2l[T3l[T4l|T5l[T6t[T7l|T8'.'T9'.'20l 
[1TJ[52] or Hartree-Fock calculations 23j as well as from tight-binding models^S, 26, 27). 
There is also an enormous amount of related work on the chemisorbed defects such as 
the hydrogen defects and chemisorbed magnetic atoms. [28| Most of the tight-binding 
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Figure 2. Splitting of the three dangling bond sp^a states of the carbon triangle, 
denoted by Vtr, and the vacancy-induced zero-mode Vtt state originating from the 
IT band. The splitting of the Vcr states is discussed in detail in Section |3] 



models have focused on the tt bands only, which is clearly inadequate due to the 
formation of the sp^a dangling bonds in the mid-gap region. The first-principles 
calculations do include all relevant states in the band structure including the sp^a 
states; however, in spite of all these works, a clear picture of the vacancy states has 
not emerged. 

In this paper, we study the electronic structure of the vacancy in graphene using 
the all-electron density functional linear augmented plane waves (LAPW) method 
along with tight-binding studies as well as the impurity Green's function (GF) 
approach to interpret the band structure. To our knowledge, this is the first all- 
electron density functional calculation for the vacancy in graphene reported in the 
literature. We have already reported the electronic structure for the mono and bilayer 
graphenes using the same method. |29j In addition to the DFT calculations, the nature 
of the vacancy-induced states is modeled from the tight-binding and Green's function 
studies, which help interpret the DFT results. 

The basic overall picture of the electronic structure that emerges from our work 
is summarized in Fig. [T] It shows the standard a and tt bands of graphene plus the 
vacancy-induced states, denoted by Vtt and Vcr, which are split due to the crystal 
field, Jahn- Teller, and the Hund's-rule interactions. The Vcr states are made out of 
the three sp^a dangling bond states, which are located on the three carbon atoms 
adjacent to the vacancy with their lobes directed towards the vacancy site. With 
their bonding partners missing, they occur in the midgap region. At the same time, 
a localized state Vtt gets introduced in the tt bands in the midgap region as well, 
the so called "zero mode" state, whose energy is exactly zero in the nearest-neighbor 
tight binding approximation. These four states, localized around the vacancy center, 
can hold eight electrons in total taking into account the spin degeneracy. The level 
structure of the vacancy-induced states is shown in Fig. [2l 

At the same time, electron counting arguments show that the vacancy releases 
four electrons to be occupied among the above localized states. These electrons include 
the three orphan sp^a electrons, one from each of the three carbon atoms adjacent 
to the vacancy, plus one orphan tt electron, whose origin may be understood in the 
following way. Focusing on the tt states now and considering a vacancy on the A 
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sublattice, the majority sublattice B has one extra atom, Nb — Nj^ — 1, so that the 
total number of tt orbitals is Nj^ + Nb, which is the same as 2iV^ + 1. Out of these, 
there is one zero-mode state and the electron-hole symmetry of the graphene lattice 
results in Na band states below E = and the same number above it. (See Fig. [7] 
for the TT band structure). So, of the 2Na -I- 1 tt electrons (one per atom), 2Na fill up 
the lower bands, leaving a lone orphan tt electron. These four orphan electrons (three 
cr and one tt) occupy the vacancy-induced states as indicated in Fig. [T] 

The remaining sections are organized as follows. In Section II, we discuss 
the results of our DFT calculations. Section III discusses the crystal-field and 
Jahn- Teller splitting of the vacancy-induced localized a states and Section IV is 
devoted to the vacancy-induced tt states. In Section IV A, we revisit the zero-mode 
theorem and in Section IV B, we present numerical results for the tt states from 
a numerical diagonalization of the tight-binding Hamiltonian before discussing the 
vacancy-induced tt states using the Green's function approach in Sections IV C and 
D. Finally the results are summarized in Section V. 

2. Density-Fiinctional Calculations 

For the density-functional calculations, we used the all-electron spin-polarized linear 
augmented plane wave (LAPW) method [50] with the general gradient approximation 
(GGA) [ST for the exchange-correlation functional. A 72-atom 6x6 supercell was used 
which included one vacancy site. The LAPW basis functions included the carbon 2s 
and 2p valence orbitals and a momentum cutoff of RK^ax — 5.2 was used, with 
approximately 3500 basis functions and about 50,000 plane waves at each k point. 
All atomic sphere radii were taken as 0.66 A. The maximum angular momentum for 
the wave function expansion inside the atomic sphere was kept at Imax = 6. Thirty 
k points in the irreducible Brillouin zone were found to be sufficient for converged 
results in the self-consistent calculations. 

Relaxed structure - First we performed a structural optimization of the lattice 
constant for pure graphene which yielded about the same lattice constant as the 
experimental value. For the vacancy calculation, the lattice constant was held fixed 
at the experimental value and a structural relaxation was performed for the entire 
structure. The optimization yielded a planar Jahn- Teller (JT) distorted carbon 
triangle around the vacancy with the carbon atoms outside the triangle relaxed by 
a much smaller amount. For the carbon triangle, we found two long bonds of length 
2.66 A each and a short bond of length 2.40 A (Fig. [3]), as compared to the 2.48 
A for the undistorted structure. In terms of the standard Jahn- Teller modes of the 
equilateral triangle, the magnitudes of the distortion are: Qo = 0.08 A (breathing 
mode), Qi = 0.166 A (symmetric bond-bending mode), and Q2 = ^ (asymmetric 
mode). [32] 

The relaxed structure for the vacancy reported in the literature varies widely. 
While some have reported planar structures, [131 HZl [TSl dni [131 E] others have found 
non-planar structures with out-of-plane displacements varying from 6z ~ 0.12 — 0.47 
A.[l2| [T9l [20I [2n [22] We find that a paramagnetic relaxation (less accurate for the 
present problem) yields a non-planar structure 5z « 0.27 A, while a spin-polarized 
calculation yields a planar structure, an observation made by Faccio et al.|17j from 
their calculations as well using the SIESTA code. We attribute this wide variation in 
the calculated relaxed structure in the literature partly to the unusual nature of the 
Vtt bound state, which falls off only as l/r, leading to a larger width of the Vtt band 
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Figure 3. The Jahn- Teller distorted planar carbon triangle obtained from the 
structural relaxation using the all-electron spin-polarized LAPW-GGA method. 

in the supercell calculations than is expected from resonance broadening due to the tt 
band continuum. 

The calculated vacancy formation energies agree much better between different 
calculations. Our result for Ey = i?(graphene -I- vacancy) — N~^{N — l)£^(graphene), 
A'^ being the number of atoms in the graphene supercell, is 7.87 eV, which compares 
well with the previous calculations [T2j [20l [15] of 7.4 — 7.8 eV as well as with the 
experimental value of 7.0 ± 0.5 eV.[55] 

Electronic structure - Fig. |3]shows the band structure, where the vacancy induced 
Vfj and Vtt states are clearly seen. The momentum points in the Brillouin zone for the 
band structure plot are defined as: K = x/\/3 + y and M = y in units of 27r3^^a^^/n 
with n = 6 for the 6x6 supercell used in the calculation. For this supercell, it can 
be easily seen by drawing both Brillouin zones that the Dirac points K and K' of 
graphene get folded into the F point of the supercell Brillouin zone, so that remnants 
of the Dirac bands are seen at the F point in Fig. |4] just above Ep (See also Fig. 
[7] for the folded graphene tight-binding tt bands for the same supercell and note the 
similarity between the tight-binding tt bands and the DFT bands, Fig. S]). Due to 
symmetry, the a and tt bands don't mix (strictly forbidden for the planar geometry, 
but also weakened significantly if the relaxed geometry is non-planar) , which leads to 
clearly identifiable vacancy-induced Vcr bands. The Vcr states originating from the 
dangling bonds are split due to the crystal field, Jahn- Teller, and exchange coupling 
as indicated in Fig. [5] and discussed in more detail in Section |31 The dispersion 
of the Vcr bands in the band structure comes from the vacancy-vacancy interactions 
in different supercells or from the k-dependent interaction with the bonding and the 
anti-bonding a bands, both effects being small. For non-planar relaxed structure, 
they should have a small resonance broadening due to the interaction with the tt band 
continuum. Three electrons occupy these states leading to the occupation Vcti fii 
Vct2 tj with the remaining fourth electron occupying the Vtt f state. 

We now turn to a description of the effect of the vacancy on the tt states. 
Basically, the vacancy introduces a sharp, resonance state Vtt in the midgap region. 
The following summarizes the discussions in Section |4l which arc important to keep in 
mind here: (i) If only NN tight-binding hoppings are kept, then the vacancy introduces 
a single localized state Vtt at i? = and of zero width called the zero-mode state, and 
its wave function decays as ^ l/r with distance in the linear-band approximation, (ii) 
Presence of the vacancy in each supercell does not affect the energy or the width of 
this state because of the result that the zero-mode wave function lives on the majority 
sublattice entirely and any changes in the minority sublattice does not affect it (in the 
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Figure 4. Spin polarized band structure of graphene with a single vacancy in a 
72-atom 6x6 supercell obtained from the density functional LAPW method. The 
vacancy induced Vct and Vtt bands are indicated in blue and red, respectively. 
Symmetry strictly forbids the admixture between a and n states for a planar 
relaxation around the vacancy, leading to flat Vcr bands (blue lines). The Vtt 
bands are not flat owing to hybridization with the continuum vr states. The Dirac 
points K and K' of the original graphene Brillouin zone get folded into the F 
point of the supercell Brillouin zone. The zero of energy is taken to be the Fermi 
energy Ep. 



supercell, all vacancies are located on the same, minority sublattice). (iii) However, 
due to the 2NN hopping as well as the exchange splitting, the energy of Vtt is different 
from zero, so that it now has a small but finite width due to resonance broadening 
with the linear tt band continuum consistent with the STM experiments. [10] (iv) In 
the supercell calculations, the Vtt state acquires an extra broadening due to the slow 
1/r decay of the Vtt wave function, because of the interaction between the supercells. 

Dirac point - In Fig. |4l the Dirac point occurs above the Ep (see the bands just 
above Ep at the F point, to which the standard Dirac points K and K' have become 
folded). For the truly isolated vacancy, the location of the Dirac point above Ep would 
mean that an infinite number of electrons are transferred from the unfilled part of the 
Dirac cones to the lone vacancy site, which is impossible. Another way of seeing this 
is to consider first an infinite lattice without the vacancy. Obviously, the Ep occurs 
at the Dirac point with zero density-of-states as usual. Now, if we introduce a single 
vacancy into the structure it can only affect the position of Ep by ~ l/iV, where N is 
the total number of atoms in the lattice, so that Ep remains unchanged for the infinite 
lattice. Of course, the electron states in the local neighborhood of the vacancy are 
modified, e. g., due to the resonance interaction with the vacancy states or due to the 
vacancy potential. The Dirac-like bands seen just above Ep at F in Fig. |4] represent 
the effect of the vacancy on the electronic structure in the local neighborhood of the 
vacancy in the supercell calculation. 
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Figure 5. Sketch of the magnetic moment for an isolated vacancy, emphasizing 
the spatial extent of the various electronic states. The Vo" electrons are highly 
localized on the carbon triangle surrounding the vacancy, while the Vtt electron 
is only "quasi-localized" with its wave function decaying only as 1/r. Hund's- 
rule exchange aligns the Vcr and Vtt electrons, producing a S = 1 state with 
the nominal magnetic moment of 2iib- This moment is however reduced by 
polarization of the tt band spins in the vicinity of the vacancy, described by an 
antiferromagnetic Kondo-likc coupling between the n bands and the localized 
Vtt and Vcr moments. The vr band polarization is about O.SfiB iii our DFT 
calculations, leading to the net magnetic moment of 1.7 fiB- 



Magnetic moment - The vacancy magnetic moment consists of two parts as shown 
schematicaUy in Fig. [5j (i) the locahzed moment coming from the vacancy states Vtt 
and Vcr and (ii) the induced moment on the band electrons in the vicinity of the 
vacancy. One can argue on general grounds that the first contribution should be 2/Lt_B 
{S = 1), while the second contribution should reduce this value somewhat due to the 
antiferromagnetic Kondo-like coupling between the localized and the itinerant band 
spins. Turning to the localized states, the vacancy leaves four electrons to be occupied 
among the Vcr dangling bond states and the Vtt zero-mode state. Of these, three 
electrons will occupy the Vcr states, so that one electron resides on each of the three 
dangling bonds of the carbon triangle. The Coulomb interaction U would prevent the 
occupation of a fourth Vcr state, so that the remaining electron is energetically favored 
to occupy the tt states. The Hund's couphng between the Vct and Vtt electrons leads 
then to a S" = 1 state with a magnetic moment of 2hb ■ This basic picture is illustrated 
in Figs. [2] and [5] and it is fully supported by the DFT bands. Fig. 21 This localized 
magnetic moment of 2fj,B is reduced due to the spin-polarization of the n bands in the 
vicinity of the vacancy. 

The spin polarization of the tt bands can occur due to two factors: (i) the 
resonance coupling with the Vtt f electron with the tt continuum bands and (ii) the 
Kondo-like antiferromagnetic interaction between the localized vacancy states and the 
continuum tt states. The first is not well described in a supercell calculation due to 
the long-range nature of the Vtt state and the second effect is intrinsically not well 
described within the band theory. 

Our DFT calculations yield a magnetic moment of about 1.7 ^b- This can be 
seen by estimating the number of holes in the small hole pocket in the two bands 
just above Ep at the T point in the spin-up bands of Fig. [4] The spin-down bands 
must contain exactly the same number of extra electrons missing from the spin-up 
bands. Without this pocket of holes, which represents the band polarization in the 
immediate neighborhood of the vacancy, the magnetic moment would be exactly 2/is, 
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corresponding to the full occupancy of Vcti ti: Vct2 t: and Vtt The existence of the 
hole pocket reduces this number. We can estimate the number n in the hole pocket 
by computing the total area of the two hole Fermi surfaces and comparing it to the 
area of the supercell Brillouin zone, which yields the value n ~ 0.15. Since the same 
number of electrons must be accommodated in the spin down bands, this would cause 
a net reduction of A^^ — iVj, by 0.30 leading to a net magnetic moment of 1.7/is. 

In the literature, the calculated magnetic moment varies widely, anywhere 
between 1.04 - IMi^b^ \M\M\MW\ES\EMMi^^^^^- Typically, 
the lower values come from calculations, where the Vtt f and Vtt J, bands overlap 
significantly. We suggest that the variation of the calculated magnetic moment in the 
literature is due to the intrinsic deficiency of the supercell method in estimating the 
TT magnetic moment due to the slow 1/r decay of the Vtt state, which produces an 
extra broadening of the Vtt state due to the supercell interaction and does not take 
into account the full anti-ferromagnetic polarization of the itinerant tt band states. 

The exchange splitting A of the Vtt state is due to its overlap with the Va states 
which are localized on the three carbon atoms adjacent to the vacancy. It may be 
estimated from the expression 

A EE E{Vtti) ~ EiVir-^) w x |*oP « 0.35 eV, (1) 

where the Hund's-rule energy is typically Jh ~ 0.9 — 1.0 eV for the atoms and 
l^oP ~ 0.4 is the combined total density of the Vtt state on the carbon triangle 
as obtained from the DFT results. The estimated exchange splitting is in agreement 
with the splitting seen in the DFT bands. Fig. 01 

Relation to the Lieb 's Theorem - The Lieb's theorem |34) states that for the 
repulsive one-band Hubbard model on a bipartite lattice and half-filled band, the 
ground state has spin S — (1/2)|A'^b — ^a|, Na {Nb) being the number of sites on 
the two sublattices. It is important to point out that the theorem holds if we consider 
only the 7r-band system and also neglect the small second-neighbor interactions that 
couples the two sublattices. Thus, with a single vacancy present, \Nb ~ Na\ = 1 so 
that according to the Lieb's Theorem we should have a net spin oi S — 1/2. However, 
in addition to the tt, we also have the a electrons. The Lieb result of 5 = 1/2 for 
the TT electrons is now coupled to the spins of the three a electrons localized near the 
vacancy, leading to the net spin S* = 1 as indicated in the summary figure. Fig. [T] We 
have already argued that the magnetic moment of 2/1^ corresponding to S* = 1 will 
be reduced due to the polarization of the band electrons in the local neighborhood of 
the vacancy. 

3. Vacancy-induced Vct states 

The essential features of the density-functional results may be understood by simple 
tight-binding considerations of the effect of the vacancy on the a and the tt bands. 
We study the a states in this section followed by the tt states in the next section. 

The description of the vacancy-induced Va states for graphene is rather simple. 
In graphene, the sp^a states are removed away from Ep due to strong interaction 
with neighbouring orbitals along the C-C bonds. However, with a vacancy present, 
the three sp^a orbitals of the three NN carbon atoms with their lobes pointed towards 
the vacancy have their usual bonding partners missing, so that they occur near Ep, 
with their on-site energies Eo- slightly below the tt orbital energies because of the s 
orbital component present in the cr states. 
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The crystal-field splitting however will lift the three-fold degeneracy. The main 
feature can be described by taking into account the 2NN hopping T between the three 
dangling bonds in the undistorted triangle, leading to the 3x3 Hamiltonian: 



diagonalization of which yields a double degenerate state &i E = T and a single 
degenerate state ai E = —2T as shown in Fig. [51 where we call this splitting the 
crystal-field splitting. The Jahn- Teller distortion of the triangle splits the double 
degenerate state further, which is described by the unequal hopping T ^ T' . 
Taking the isosceles-triangle relaxation found in our DFT results, two of the three 
hopping terms are modified into T' as indicated in Fig. |6l From the DFT band 
structure, we find that T w 1.6 eV, while T' « 1.2 eV. The new eigenstates are: 
£'cti,o-2 = 2~^(-T=F V8T'2 _|_ j^2) and E^^ = T with the corresponding (unnormalized) 
wave functions ^-1^2 = ((-T ± V8T'2 _^ y2) /T', 1, 1) and *3 = (0,-1,1). This simple 
model suggests a Jahn- Teller distortion of the carbon triangle surrounding the vacancy. 

The Jahn- Teller interaction is of the type E ® e (both electronic and vibrational 
states are doubly degenerate) in a trigonal (D3/1) symmetry. With this lattice 
distortion, the trigonal symmetry is broken. The distortion removes the double 
degeneracy and the two states (shown in Fig. [2] as Va2 and Va^) are now split 
by the amount 2^1 (3r - VST'^ -|- T^) w 4(T - r')/3 w 0.55 eV. Since there are only 
three electrons available to the Vct states, Vai is occupied with two electrons, while 
the lone remaining electron occupies the Va2 state. The spin-degeneracy is removed 
by the Hund's coupling with the electron occupying the Vtt state, producing the spin 
structure indicated in Fig. [21 The wave function ^'2 corresponding to the Va2 state 
shows that the maximum weight (~ 66%) comes from the sp^cr dangling orbital of the 
apical atom of the carbon triangle, which is consistent with the spin density plotted 
in Fig. [HI The Jahn- Teller distortion is actually dynamic, with the carbon triangle 
tunneling between three equivalent minima on the adiabatic potential surface of the 
E ® e Jahn- Teller problem, an issue we discuss elsewhere. [35] 

4. Vacancy-induced Vtt states 

In this Section, we discuss the origin of the localized state - the so-called "zero- 
mode" state - near the band center of the tt bands. Understanding of the origin 
and the "quasi-localized" nature of the zero-mode state is an essential part of the 
interpretation of the full band calculation using the density-functional theory. 

If only the NN interactions are present, the "zero-mode" state is a localized state 
with energy exactly at the band center. If in addition the higher-neighbour interactions 
are also present but not too large, as is the case for graphene, f29^ then the localized 
state turns into a sharply-peaked resonance state owing to its overlap with the tt bands 
and occurs not too far from the band center. 

^.1. The existence of the zero mode state 

According to the zero- mode theorem. [25j which is in fact valid for any bipartite lattice 
with NN interactions, whenever there is an imbalance in the number of atoms in the 
two sublattices in a bipartite lattice, viz., n = Nb — Na > 0, there are n number 




(2) 
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Figure 6. Spin density — at different sites in graphene witfi a vacancy as 
obtained from the density-functional calculations. Green (blue) circles indicate 
positive (negative) values, with the area of the circle being proportional to the spin 
density. The spin moments on the carbon atoms other than the vacancy triangle 
are due to the 7r electrons, which are long-ranged due to the slow 1/r decay of 
the Vtt state. The hopping integrals T and T' between the sp^a orbitals on the 
carbon triangle adjacent to the vacancy has reference to the model discussed in 
Section [3] 



of degenerate solutions with the eigenvalue es (the on-site energy of the majority 
sublattice), with the wave functions residing entirely on this sublattice. This can be 
seen from the following simple considerations as an alternative to Pereira et al.'s proof 
which used the rank- nullity theorem in linear algebra. [25] 

We begin with the conjecture that there are some solutions where the wave 
functions live completely on the majority sublattice (B) and proceed to find them. 
Thus we have 

where is a vector in the B sublattice of dimension Nb and there is null contribution 
from the A sublattice to the wave function. It will be clear from the following 
discussion that for the theorem to hold, the B sublattice Nb x Nb Hamiltonian is 
restricted to the diagonal form 

Hbb = es/, (4) 

and there are no restrictions on the remaining parts of the Hamiltonian. The specific 
form oIT-Lbb means that there is no site disorder, nor is there any interactions between 
the atoms on the B sublattice (hence it will fail if interactions beyond the NN are 
present, which will produce a non-diagonal Hbb)- However, such restrictions need 
not apply to the A sublattice, so that the Na x Na Hamiltonian Haa for the minority 
sublattice can have diagonal disorder and also there is no restriction on the form of 
T-l BA as well. This means that the A sublattice atoms can interact between themselves 
and with the B sublattice atoms as well without invalidating the theorem. 
The wave function "i! b thus satisfies 

= E'fB (5) 

^L*B = 0. (6) 

The first of these equations tells us that if the conjectured solutions of the form 
(^"3, 0) exist, then they must have the enegy E — €b and there would be at most Nb 
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number of such degenerate solutions; the equation does not place any constraints on 
the individual components of 

Turning to Eq. [51 there are Ng components of g to be determined but only 
Na < Nb equations are there to determine them. This means that the solutions 
can not be fully determined. However if we specify Nb — Na components of 
then the remaining components are uniquely determined as linear functions of these 
components. These solutions are therefore of the form 

*B = {<Pl,<p2, --Anb-Na^ f2, Ina)! (7) 

where we can choose the 0,;'s arbitrarily and the /^'s are then just linear combinations 
of 4>iS {fi = X^j^i ^'^'^ij'Pj^ where the expansion coefficients are determined by 

"H^BA i'^ -f^l- ©■ Thus the number of linearly independent solutions is given by 
the number of ways we can choose linearly independent solutions in the subspace 
(01:02, ■■■,4'Nb-Na)j which is clearly Nb — Na- This proves the conjecture and the 
theorem. 

It is easy to see why the theorem is not valid if there is on-site disorder on the 
majority sublattice or interactions beyond the NN, which introduces off diagonal terms 
in Hbb- So, Eq. |4]is not true anymore. This means that Eq. [5]puts constraints on the 
components of V&s in order to satisfy the eigenvalue problem and as a result Eqs. [5] and 
[6]can not both be satisfied simultaneously. For example, if we use the form Eq. [7] which 
satisfies Eq. [HI we are only left with the freedom to choose 01, 02, (^Nb-Na ^^'^ this 
is not enough to satisfy the eigenvalue problem of Eq. [5l There is no such problem 
if Hbis — EbI, since any vector {(pi,(p2, ■■■,<Pnb-Na'j fi^ f2i fNA) is a solution with 
E = eB. 

The theorem has an important bearing on the results of the supercell calculations, 
both tight-binding and density functional. In these calculations, the vacancies are 
repeated in each supercell, connected by the superlattice translational vectors, and 
hence are all located on the same sublattice, which forms the minority sublattice. If 
n is the number of supercells in the crystal, then this is also the imbalance in the 
number of atoms in the two sublattices n — Nb — Na- According to the theorem, 
there should be n zero-modes in the Brillouin zone, which is also precisely the number 
of Bloch momentum points in the Brillouin zone. These states thus show up in the 
form of a dispersionless band in the tight-binding supercell calculations a.t E = 0. 

If hopping beyond the NN is present or if the on-site energies of the different 
atoms are different, then the theorem does not hold. However, the hopping beyond 
the NN in graphene is small and the on-site energies are negligibly different on sites 
close to the vacancy as obtained from the DFT calculations. Because these effects are 
small, a clearly identifiable, nearly-dispersionless zero-mode band is found in the DFT 
calculations as seen from Fig. |4]as well as in the higher-neighbour tight-binding results 
[Fig. [7], where the zero- mode band is marked by the red dots. 

4-2. Tight-binding results: Direct diagonalization of the Hamiltonian 

In order to further understand the formation of the zero-mode states, we have 
studied the vacancy tt bands with the standard tight-binding model of the Pz 
orbitals containing both the nearest neighbour (NN) and the second- neighbour (2NN) 
interactions. In particular, we look for the vacancy-induced zero-mode states discussed 
in the previous Subsection. 
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Figure 7. Tight-binding band structure obtained from the Hamiltonian Eq. |8] 
for the 72-atom 6x6 supercell both with and without a vacancy. The band 
structure without the vacancy is shown in the left paneh The Dirac points 
K and K' of the original graphene Brillouin zone get folded into the F point 
of the supercell Brillouin zone and the two linear Dirac bands are clearly seen 
in the left panel. The middle and the right panels show the zero-mode states 
(red dots) with and without the second neighbour interaction t' . In the NN 
tight-binding approximation (middle panel), all zero-modes have the same energy 
and live exclusively on the majority sublattice, while with the second neighbor 
interaction, the zero-mode states do have a band dispersion and leak into the 
minority sublattice as well. The sublattice with the vacancy atoms is labelled A 
and the total number of states in different bands (not counting spin degeneracy) 
is shown on the right, with N,^ and A'^^ denoting the total number of atoms in 
the two sublattices in the entire crystal. 

The tight-binding Hamiltonian is 

nTB = -tJ24a^J-+t' E C,W+i/.C., (8) 

where —t and t' are the NN and the 2NN interactions with the signs chosen such that 
t,t' >0 {tK 2.91 eV and t' « 0.16 eV for graphene [13] )• 

The band structures and the densities-of-states are shown in Figs. [7]and|Sl The 
electron counting in the band structure Fig. [7] is as follows. Both the lower and the 
upper bands contain in total (integrated over the Brillouin zone) Na states each, while 
the zero-mode band contains Nb — Na states, making a total of Na + Nb states, as 
it must be the case. We have one tt electron per site present in the system, so that 
taking spin into account, the entire lower subband is full, while the zero-mode states 
are half full. For the single vacancy {Nb — Na = 1), this leads to a single occupied 
electron in the zero-mode states, resulting in 5 = 1/2 in agreement with the Lieb's 
Theorem f34|. 

As discussed in the previous Section, if only the NN interactions are present, we 
should have Nb — Na number of zero-mode states at i? = exactly. This is why 
the zero-mode band in the middle panel of Fig. [7] is completely flat. However, if the 
2NN interactions are also present, then the energies of the zero- mode states are not 
guaranteed to be the same and we see a spread in the energy of these states, which 
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Figure 8. The tight-binding n densities of states of graphene with a vacancy 
with NN interacations only (left) and with both the NN and the 2NN interactions 
present (right) as obtained from the tight-binding Hamiltonian Eq. [S] 



shows up as a dispersion in the zero-mode band, as seen in the right panel of Fig. [71 
Here, the vacancy site was modeled by simply removing a lattice site, 
corresponding to the vacancy potential Uq — oo. In a real material, however, Uq 
is large but finite. The effect of a finite Uq is that (a) It causes the zero-mode state 
to occur shghtly below the mid-gap {E = 0) and (b) The sharp zero-mode state turns 
into a resonance state due to interaction with the continuum tt bands. This is best 
described with the Green's function approach discussed in the next Subsection. 

4-3. Impurity Green's Function and the zero-mode state in the tt bands 

In this Section, we study the effect of a single impurity on the n electron states by 
studying the Dyson's equation and show the emergence of the zero-mode state as the 
strength of the impurity potential Uq is gradually increased. For the vacancy, this 
potential is large but finite, so that the results obtained in this Section are helpful in 
understanding the nature of the zero-mode state in the actual structure with a finite 
vacancy potential. 

The wave function of the zero-mode state is obtained from the Lippmann- 
Schwinger equation. Consistent with the previous results. [2^ I36j we find that (a) 
The zero-mode state consists of wave functions from the majority sublattice only and 
(b) It is quasi- localized decaying as 1/r as a function of distance from the vacancy in 
the limit of the linear-band approximation. In addition to these known results, our 
analysis allows us to (a) obtain the oscillatory phase factors in the zero-mode wave 
function due to the interference of the two Dirac points and (b) compare the linear- 
dispersion results with the full tight-binding band result by computing the GFs for 
large distances in both cases. 

The vacancy is modeled by adding an on-site perturbation V to the unperturbed 
NN (NN) tight-binding (TB) Hamiltonian, so that 

n = no + v, (9) 

where Hq — ^t^cl^Cjp + H.c, ia being the site-sublattice index, and 

V = Uocl^coA, (10) 
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Figure 9. On-site GF Gq^ oa(^) ^'-"^ bands computed from the Horiguchi 

method and the energy of the resonance state, indicated by the black dot, obtained 
from the Dyson's equation: UoFo{E) = 1. As Uo — > oo, the solution moves to 
i? — > leading to the sharply-localized zero-mode state at the band center. 
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Figure 10. Change in the total sublattice DOS due to the addition of the 
impurity with Uo = 5t as computed from the factors multiplying the coefficient 
"1/Af" in Eqs. 1171 and 1181 The emergence of the zero-mode state on the B 
sublattice at i? = is clearly seen, which grows into a 5 function as Uo — > oo. 
The remaining changes in the DOS go to zero as 1/A'^, except for the (unimportant) 
bound state on the A sublattice (occurring at E/t 6 in the top figure), which 
becomes a iS-function bound state with i? — > oo in the limit of Uo c>d. 



where Uq is the strength of the potential due to the impurity on the A-sublattice in 
the central cell. The vacancy corresponds to the value of C/q — ^ oo. 

The key quantity of interest here is the unperturbed GF, G^{E) = {E+irj—T-Lo)^^ , 
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the calculation of which we have discussed in our earlier paper where we studied the 
RKKY interaction in graphene. [33 usual, the imaginary part of the GF contains 
the information about the density-of-states: p°{E) = -Tr-^Im G°{E). The GF G{E) 
in the presence of the perturbation will be obtained from the Dyson's equation. 

Since we will be interested in the local density-of-states (LDOS) on the various 
carbon sites and how they are modified by the presence of the impurity, we will need 
to calculate the real-space matrix elements G^^jp{E) = {ia\G^ {E)\j /3) . This may 
be done by going to the momentum space and defining the Bloch functions for the 
electrons \ka) — A^~^/^ e*'^''"'" with fia — Ri + Ta being the position vector 
of the a-th atom in the i-th unit cell. The unperturbed Hamiltonian "Ho in this 



basis set becomes "H,- = ( '^'"'^^ ] , where f{k) -t [e'^-^^ + e'^''^^ + e'^''^^ 

/*(fc) 



and di, d2, and dz are the positions of the three nearest neighbors. Diagonahzation 
of the Hamiltonian yields the eigenenergies E{k) — ±|/(fc)|, which when expanded 
around the Dirac points lead to the usual linear band structure E{q) — ±Wi?|(j|, where 
q = k — Kd is the deviation from the Dirac point in the Brillouin zone. Here the 
Fermi velocity vp — ita/2, where 'a' is the carbon-carbon bond length. Note that 
unlike our previous work. [37] vp here is defined to be a positive quantity, since 't' is 
positive. 

The real-space GFs are conveniently obtained by first calculating the momentum- 
space GF, which can be easily shown to be G^{kE) = (fca|G°(£')|fc/3) = {E + irj + 
H^)((i?+i77)^ — |/(fc)P)^^. A Fourier transform then yields the real-space unperturbed 
GF, viz.. 



GL,,p{E) = ^Y.^~^'^''"~'"'^G%{kE), (11) 

k 

which can be calculated by simply a brute-force summation over the Brillouin 
zone. It can also be computed by a second method using the Horiguchi recursive 
technique. [3H1 133 However, the latter, although computationally fast, has convergence 
problems [39] for distances \Ri ~ Rj\ > la or so, so that this is a better method to use 
only for smaller distances. 

The perturbed GF is related to the unperturbed GF through the Dyson's 
equation: G = + G'^VG. Using the localized form of the impurity potential, 
Eq. [TUl and taking the matrix elements, we immediately get for the real-space GF, 
the result 

G.„,,^(i?) = G\^^^f,{E) + C/o X G\^^^j,{E)G^A,,fi{E). (12) 
We are specifically interested in the on-site GFs with a — fi and Ri = Rj , which 
give the LDOS on the A and B sites at distance ria = Ri -\-Ta from the impurity site. 
Eq. (|T2|) is easily inverted to yield the perturbed G{E) in terms of the unperturbed 
G^iE), viz.. 

The LDOS at different sites may be obtained by taking the imaginary parts of 
the diagonal elements of the GF: pia{E) — — Tr^^Im Gia,ia{E). It immediately follows 
from Eq. |T3]that the LDOS at the impurity site has the much simpler form 

''^^''^ {l-UME)f + {.U,p,{E)r ^''^ 



Electronic structure of the substitutional vacancy in graphene: Density-functional and Green's function studieslG 



where po{E) = — Tr^^Im Gq^q^(£') is the unperturbed LDOS at the central site, 
which is of course the same for every site in unperturbed graphene, and Fq{E) = 
Re Gq^ oa(^)- Note that we have defined here po{E) to be the unperturbed density- 
of-state per sublattice per spin (which is independent of the sublattice or the cell 
index) and pia{E) is the corresponding perturbed quantity for the ia site. 
The resonance condition follows from Eq. I14| viz., 

l-f/oFo(£;o) =0, (15) 
where Eq is the resonance energy. The graphical solution of this equation is shown in 
Fig. [9l There are four solutions for Eq: The two solutions at Eq — ±t do not produce 
much change in the DOS, as may be inferred from Eq. [T31 due to the diverging 
density-of-states po{E) there, and the bound state with Eq — > Un for large Un is 
inconsequential because it is removed to oo, which then leaves the sole resonance state 
indicated by the black dot in Fig. [9l Its energy goes to zero in the limit Uq ^ oo and 
the oscillator strength to one, producing the zero-mode state for the vacancy. 

The total DOS in the presence of the perturbation may be computed by taking the 
trace of Eq. for the entire lattice. Using the identity J^ia oa(^)^oa lai^) = 
— cJGq^ Qj^{E)/dE and some tedious algebra, the result is 



. , 1 ^^ -Uo[Uopo{E)Fl,{E)+p'^{E){l-UME))] 
Ptot(i^) = ^P.{E) + - X ^^_uME)f + i.U,p,i^E)Y • ^''^ 

Similarly, by taking the trace of Eq. ([T^ over the cell index only, the individual 
sublattice DOS may be found, which for the A sublattice reads as 

Pa(E) = po(E) + — X — ——^ , 17) 

N t:[{1 - UoFo[E))'' + [7rUoPo[E))'^\ 

where £,{E) = {i/N)'^i^[G'^^{kE)]'^ and the densities of states are, again, per 

sublattice and per spin. A similar expression for pb{E) reads 

Pb{E) = po{E) + 1 



-Uo[{l - UoF^{E)){ttp'^{E) - Imi{E)) + ttU^p^{E){F{,{E) + Rei{E))] 



.(18) 



^[(1 - UoFo{E)f + (irUoPoiE))^] 
It can be verified from Eqs. - that Piot{E) = Pa{E) + Pb{E), so that these 
equations are consistent. 

The numerical results are summarized in Figs. [TUl [TTl and [T^l The factors 
multiplying the 1/N in Eqs. [TTl and \TE\ are the changes in the DOS Apa{E) and 
Apb{E) introduced by the impurity potential, which are shown in Fig. [TUl Fig. [TH 
shows the emergence of the zero-mode in the density-of-states with E — and that 
this state resides entirely on the B sublattice in the limit C/q — )■ oo. Fig. [121 shows the 
LDOS on the impurity site (poa) and on the nearest (pob) and the next nearest sites 
(Pia)- 

The width of the resonance peak increases with the resonance energy Eq . Keeping 
the linear term in the expansion of Fo{E), viz., Fo{E) « Uq^ + Fq{Eq){E — Eq), Eq. 
P6I) yields the Lorenzian lineshape 

PtotiE)^2pQiE) + J^^^-^___, (19) 

where the resonance width is given by 

r = -npoiEo)/F^iEo). (20) 
The width is zero if Eq = Q and increases with energy as shown in Fig. (fT3|) . 
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Figure 11. Total density-of-states for the A sublattice pa{E) (top) and 
the B sublattice pg{E) (bottom) for different values of the impurity potential 
Uo/t = 0,2, and 5, denoted by the black dashed, black solid, and red dashed lines 
respectively. These results are obtained from Eqs. ll7l and[T8lbv using Af = 20 for 
the purpose of plotting. The figure shows the evolution of the zero-mode state at 
E = 0, which lives completely on the B sublattice in the limit of Uo oo, i.e., 
opposite to the sublattice in which the vacancy is introduced. 



4- 4- Impurity state wave function 

In this section, we study the impurity state Vtt wave function from the Lippmann- 
Schwinger equation. The analysis allows us to obtain the well known 1/r decay of 
the vacancy state; however, in addition we also obtain the oscillatory behavior of the 
wave function due to the interference between the two Dirac cones. 

Our starting point is the Lippmann-Schwinger equation j^I') ~ I^J/") + G'^V\'i/), 
where j^*^) is the unperturbed state. For the localized impurity potential on the 
central A site, V — Uo\OA){OA\, the Lippmann-Schwinger equation leads to the wave 
function 

*,„M.oW-< + ^%^. (21) 

We are interested in the low-energy behavior, since that's the region where the 
resonance state gets introduced by the impurity as seen from Fig. [Q] The GFs for the 
full tight-binding band structure as well as for the linear bands were computed in our 
previous work.[33 For the linear bands and in the low energy limit, the results are 

ZTTVp VF 
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Figure 12. Local density-of-states at the impurity site p\A (top), the NN site 
Pos (middle), and the next NN site p\A (bottom) obtained from Eqs. 1131 and 
1141 for different strengths of the impurity potential C/oA = 0,2, and 5, denoted 
by black dashed, black solid, and red dashed lines respectively. As (7o — ^ oo, the 
top LDOS goes to zero (except for the bound state beyond the top of the band 
whose energy goes to oo), and the zero-mode state lives only on the B sublattice, 
as indicated from the middle and the bottom panels. The prominent zero-mode 
peak in the middle panel for (7o /t = 5 will develop into a 5- function peak at = 
as the impurity potential C/q ~^ oo- 



'2Trvj, 



-iEr ^ 

VF 



(22) 



where Ac is the unit ceh area in graphene, Kq and Ki are the modified Bessel functions 
of the second kind and r is the distance vector between the two atoms: r = riA — "Tqa 
for the first GF and r = ris — r^A for the second. The multiphcative factors are 
/3 = e*^ *" -|- e*^ which is a real number for the graphene lattice and 

a = e— /3(e»(J?---9.) _ ^^(K' -r+e^)^^ (23) 
which is purely imaginary and the polar angle 9r = tan~^(j//x) is defined with the x 
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Figure 13. Resonance width of tlie zero-mode state (Vtt) as a function of the 
resonance energy Eq. Both Eq and F are in units of the NN hopping, with the 
value t ^ 2.56 eV |29| if only the NN hopping is kept. 



direction taken to be along the vector K' — K, which are two adjacent Dirac points 
in the Brihouin zone. 

Using the small z expansion for the Bessel functions: [40] Ko{z) — — ln(z/2)— 7 and 
Ki{z) = + 2"izln(z/2) + (7 - l/2)z/2, where 7 ss 0.577 is the Euler-Mascheroni 
constant, we find the Bessel functions entering the expressions for the GFs (Eq. [22| 
to become, in the low energy limit: 

, iEr ^ iiT , , 2vp 
i^o(- — ) = yS.gn(i^)+ln^-7, (24) 

iEr vp iEr iEr I iEr 

Ki{ ) = --TT-T^ — ln(-:^ — ) - (7 - 7;)7; — • 

Vp lEr Zvp 2vp 2 2vp 

Plugging these into Eq. [22l taking the E ~> limit , and retaining the lowest-order 
terms in energy, we find the results 

^ ^4Jma VP E^r 9-1 ,> l^kx 

- S;;^^^Sni?]. (25) 

There is no guarantee that these results, calculated for the linear bands extending 
to infinite energy, should agree even for low energies with the GFs for the actual 
bands, e.g., as obtained with the tight-binding band structure. Certain elements must 
exactly agree at low energies, for example, the imaginary part of Gq^ q^, which yields 
the density-of-states, since at low energies, it is controlled by the Fermi velocity vp 
alone. We nevertheless find that the expressions Eq. [25] do agree quite well with 
the tight-binding GFs, the agreement becoming better with increasing distances. A 
comparison between the low-energy GFs Eq. [25] and the full GFs for a typical case is 
shown in Fig. [TJ] which also illustrates the symmetry of the GFs. A notable exception 
is the real part of the on-site GF Gg^ oa(-^)> where the substitution of r = in Eq. [25] 
yields a divergent result. However, we find that the tight-binding GF in this case can 
be fitted to the expression Eq. [25]for Gq^ Q^(i?), provided we use the value r w 0.6 a 
instead of r = 0. 

We note that the symmetry properties of the above GFs Eq. [51] are consistent 
with the results [38] that follow from the particle- hole symmetry and valid for all 
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Figure 14. Illustration of the symmetry of the GF and its low energy behavior, 
calculated using the Horiguchi method and the tight-binding band structure. 
Dashed and full lines denote the real and the imaginary parts, respectively. For 
the upper figure, the distance vector of the atom with respect to the impurity 
is given by r = r^^ — rpA = 2\/3a{l,0), while for the lower figure, it is 
r = fig — tqa = 2a(0, 1), where the coordinates are indicated in Fig. |6l The 
points near E = are the low-energy results for the linear bands as given by Eq. 
[25] 



energies, viz., Re G^^^^iE) = t Re G%^,^{E), Im G^^^^iE) = ± Im G%^^JE), and 
^ia jp^^) ~ ^% ia(^)^ where the upper (lower) sign is for a — The symmetry 

properties are illustrated for specific cases in Fig. [TH 

The nature of the impurity state immediately follows from the Lippmann- 
Schwinger expression Eq. 1211 First of all, notice an important point from the 
expression for the GF Eq. [25l viz., that all GFs vanish at £' = except for the 
real part of G'^g q^, which is finite and decays as 1/r. This is precisely what leads to 
the property that the zero-mode state resides on the majority sublattice B only and 
its wave function decays inversely with distance. These features are true if only the 
NN interactions are present on the graphene lattice. If second NN interactions are 
present, then there is no electron- hole symmetry and the behavior of the GFs near 
the resonance energy differs from Eq. [25l The form of the GFs for the latter case is 
such that both sublattices contribute to the resonance state near E ^ 0, an issue that 
is discussed in detail elsewhere. [iT] 

Returning to the Lippmann-Schwinger equation Eq. [5T] and inserting into it the 
low-energy expansion for the GFs (Eq. [25)1 and then taking the limit of the resonance 
energy i?o = 0, it can be easily seen that as E'o ^ in the limit Uq ^ oo, the impurity 
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Figure 15. Square amplitude I'tp of tiie zero-mode state on tiie B sublattice 
along the zigzag and the armchair directions computed from the Lippmann- 
Schwinger result Eq. 1271 using Green's Functions obtained for the (a) full TB 
bands (black solid lines) and (b) linear bands (analytical expression, Eq. I28I I 
(red dashed lines). Circles indicate the same quantity computed from the direct 
diagonalization of the tight-binding Hamiltonian on a finite lattice consisting of 
a single vacancy in a 512-atom supercell. 



wave function follows the behavior 




This is an important result, which states that in the NN approximation, only the 
B sublattice component survives for the zero-mode state, it being the stronger infinity. 
The surviving component is found to be simply proportional to the real part of the 
inter-sublattice GF, 

(X Re G%^aA{Ea ^ 0), (27) 

since its imaginary part vanishes. Using Eq. [2S]and evaluating Im a from Eq. [231 we 
finally get the desired result 

N ^ ^ - - 

= — sin[(if - K') ■ f/2 - 0r] cosHK + K') ■ r/2 - 7r/3], (28) 
r 

where we have suppressed the cell index i, is a constant, r is again the actual 

distance vector of the B site with respect to the impurity position, and the two 

Dirac points in the Brillouin zone may be taken as: K — 27ra^^3-'^/^(— 1, V3) and 
i^' = 2^-13-3/2(1, 73). 
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Eq. [21] is the central result of this Subsection that describes the l/r decay of 
the vacancy-induced Vtt state along with the phase factors. The long-range nature 
l/r of the wave function is well-known. [55] but the oscillatory factor due to the 
interference effect of the two Dirac points is new. The same kind of interference is also 
present in the oscillations of the RKKY interactions. [42] [37] The wave function is not 
square integrable because we used the linear band structure, but it will be if we take 
the full band structure into account. Eq. [28] nevertheless describes the gross features 
of the zero-mode state. The wave function changes sign along different directions, e. 
g., it changes sign along the zigzag direction but not along the armchair direction. The 
kinetic energy gained by the delocalization of the wave function is exactly cancelled by 
anti-bonding components present in the wave function, so that its energy still equals 
the on-site energy in spite of the delocalization. The calculated wave function for the 
zero-mode state is shown in Fig. [TS] We note that a recent study has shown that the 
l/r decay of the vacancy state remains unchanged even when a repulsive Coulomb 
interaction is included in the tight-binding Hamiltonian. [15] 

5. Summary 

In summary, we have studied the electronic structure of graphene with a single 
substitutional vacancy from density-functional calculations using the all-electron 
LAPW method and interpreted the results with the help of the tight-binding model 
and the impurity Green's Function approach. We find that the vacancy induces four 
localized states, viz., three Vcr dangling bond states on the carbon triangle and one 
Vtt resonance state. The dangling bond states cause a Jahn- Teller distortion, which 
we found to be a planar distortion of the carbon triangle. Hund's coupling between 
these electrons would then produce the S — 1 state at the vacancy center as indicated 
in the summary figure Fig. [TJ The magnetic moment has two components: (i) The 
component coming from the localized vacancy states Vcr and Vtt and (ii) An 
opposite component of several tenths of /xb coming from the spin-polarization of the 
continuum tt band states in the vicinity of the vacancy. The second part is not 
well described in the supercell band calculations due to the slow l/r decay of the 
"quasi-localized" Vtt wave function. This long-range decay also means that in an 
experimental sample it is only for the extremely low vacancy concentration that the 
truly isolated vacancy limit is reached and as a result the magnetic moment is expected 
to be dependent on the vacancy concentration. 

In addition to the density-functional calculations, we also studied the formation 
of the Vtt state in detail from the impurity Green's function approach for the isolated 
vacancy, which provided important insight in the interpretation of the results of the 
band calculations and the formation of the zero-mode states in the tt bands. This 
zero mode state is a slowly-decaying localized state that lives mostly on the majority 
sublattice. It spreads into the minority sublattice (the one containing the vacancy) 
and becomes a resonance state due to the second and the higher-neighbor interactions 
as well as the finite strength of the vacancy potential. The Green's function approach 
provided a sinusoidal phase factor associated with the Vtt wave function described 
by Eq. [28] In addition to the understanding of the vacancy electronic structure, 
our work provides important insight necessary for the understanding of impurities in 
general such as iron and cobalt dopants and other complex defects. 



Electronic structure of the substitutional vacancy in graphene: Density-functional and Green's function studies23 



Acknowledgments 

This work was supported by the U. S. Department of Energy through Grant No. 
DOE-FG02-00ER45818. 

* Permanent Address: Department of Physics, Indian Institute of Technology 
Madras, Chennai 600036, India 

** Permanent Address: Institute of Nuclear Sciences, Vinca, University of 
Belgrade, P. O. Box 522, 11001 Belgrade, Serbia 

References 

[1] Castro Neto A H, Guinea F, Peres N M R, Novoselov K S and Gcim A K 2009 Rev. Mod. Phys. 
81 109 

[2] Abergel DSL, Apalkov V, Berashevich J, Ziegler K and Chakraborty T 2010 Adv. Phys. 59: 
4 261 

[3] Bhattacharya A, Bhattacharya S, Majumder C and Das G P 2010 J. Phys. Chem. C 114 10297 
[4] Hashimoto A, Suenaga K, Gloter A, Urita K and lijima S 2004 Nature (London) 430 870 
[5] Ewels C P, Telling R H, El-barbary A A, Heggie M I and Briddon P R 2003 Phys. Rev. Lett. 
91 25505 

[6] Nortlund K, Keinonen J and Matilla T 1996 Phys. Rev. Lett. 77 699 

[7] Esquinazi P, Spemann D, Hohne R, Setzer A, Han K -H and Butz T 2003 Phys. Rev. Lett. 91 
227201 

[8] Kushmeric J G, Kelly K F, Rust H-P, Halas N J and Weiss P S 1999 J. Phys. Chem. B 103 
1619 

[9] Barzola-Quiquia J, Esquinazi P, Rothermel M, Spemann D, Butz T and Garcia N 2007 Phys. 
Rev. B 76 161403(R) 

[10] Ugeda M M, Brihuega I, Guinea F and Gomez-Rodn'guez J. M 2010 Phys. Rev. Lett. 104 096804 
[11] Cervenka J, Katsnelson M I and J. Flipse C. F 2009 Nature Phys. 5 840 

[12] El-Barbary A A, Telling R H, Ewels C P, Heggie M I and Briddon P R 2003 Phys. Rev. B 68 
144107 

[13] Yazyev O V and Helm L 2007 Phys. Rev. B 75 125408 

[14] Choi S, Jeong B W, Kim S and Kim G 2008 J. Phys.: Cond. Mat. 20 235220 
[15] Singh R and Kroll P 2009 J. Phys.: Cond. Mat. 21 196002 

[16] Yang X, Xia H, Qin X, Li W, Dai Y, Liu X, Zhao M, Xia Y, Yan S and Wang B 2009 Carbon 
47 1399 

[17] Faccio R, Fernandez- Werner L, Pardo H, Goyenola C, Ventura O N and Mombru A W 2010 J. 

Phys. Chem. C 114 18961 
[18] Telling R H, Ewels C P, El-Barbary A A and Heggie M I 2003 Nature Mater. 2 333 
[19] Lehtinen P O, Foster A S, Ma Y, Krasheninnikov A V and Nieminen R M 2004 Phys. Rev. Lett. 

93 187202 

[20] Ma Y, Lehtinen P O, Foster A S and Nieminen R M, 2004 New J. Phys. 6 68 
[21] Lim D-H, Negreira A S and Wilcox J 2011 J. Phys. Chem. C 115 8961 

[22] Dai X Q, Zhao J H, Xie M H, Tang Y N, Li Y H and Zhao B 2011 Euro. Phys. J. B bf 80 343 
[23] Forte G, Grassi A, Lombardo G M, La Magna A, Angilella G G N , Pucci R and Vilardi R 2008 

Phys. Letts. A 372 6168 

[24] Palacios J J and Yndurain F, |arXiv: 1203,_64 "85l vl 

[25] Pereira V M, Lopes dos Santos J M B and Castro Neto A H 2008 Phys. Rev. B 77 115109 
[26] Hjort M and Stafstrom S 2000 Phys. Rev. B 61 14089 

[27] Palacios J J, Fernandez-Rossier J and Brey L 2008 Phys. Rev. B 77 195428 

[28] See, for example, Ref. [13; and Wu M, Cao C, and Jiang J Z 2010 New J. Phys. 12 063020 

[29] Nanda B R K and Satpathy S 2009 Phys. Rev. B 80 164530 

[30] Blaha P et al, WIEN2k, ^' An Augmented Plane Wave + Local Orbitals Program for Calculating 
Crystal Properties" (Karlheinz Schwarz, Techn. Universitat Wien, Austria, 2001) ISBN 3- 
9501031-1-2 

[31] Perdew J P, Burke S and Ernzerhof M 1996 Phys. Rev. Let. 77 3865 

[32] Grosso G and Parravicini P 2000 Solid State Physics (Acedemic Press, London). 

[33] Thrower P A and Mayer R M 1978 Phys. Stat. Sol. A 47 11 

[34] Lieb E H 1989 Phys. Rev. Lett. 62 1201 

[35] Z. S. Popovic, B. R. K. Nanda, and S. Satpathy, In Preparation. 



Electronic structure of the substitutional vacancy in graphene: Density-functional and Green's function studies2A 

[36] Pereira V M, Guinea F, Lopes dos Santos J M B, Peres N M R and Castro Neto A H 2006 Phys. 
Rev. Lett. 96 036801 

[37] Sherafati M and Satpathy S 2011 Phys. Rev. B 83 165425;2011 Phys. Rev. B 84 125 416 

[38] Horiguchi T 1972 J. Math. Phys. 13 1411 

[39] Berciu M 2009 J. Phys. A: Math. Theor. 42 395 207 

[40] Gradshtcyn I S and Ryzhik I M 1980 Tables of Integrals, Series, and Products (Academic Press, 

New York), Sec. 8.446 
[41] Sherafati M and Satpathy S 2011 Phys. Stat. Sol. B 248, 2056 
[42] Saremi S, 2007 Phys. Rev. B 76 184430 
[43] Chang Y C and Haas S 2011 Phys. Rev. B 83 085 406 



